wfc3_dash on DASH data¶By the end of this tutorial, you will:
wfc3_dash.Introduction
1. Imports
2. Download relevant data
3. Run DASH
4. Plot original IMA vs. DASH pipeline science result
5. Conclusions
Additional Resources
About the Notebook
Citations
This notebook is the first in a new Drift And SHift (DASH) pipeline workflow developed to ease the process of reducing DASH data. The pipeline is customizable, able to be changed according to scientific goals of the user, and this first tutorial walks the user from data download to a finished product ready for science analysis. For more information, see Momcheva et. al 2016 and WFC3 ISR 2021-01.
This notebook assumes you have created the virtual environment in WFC3 Library's installation instructions.
We import:
%matplotlib notebook
import os
from glob import glob
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
from astropy.io import ascii
from astroquery.mast import Observations
from drizzlepac import astrodrizzle
from reduce_dash import DashData
obsTable = Observations.query_criteria(proposal_id=['15238'], obs_id=['IDNM0J030'])
Get the full list of products associated to the table and restrict the list to IMA files.
product_list = Observations.get_product_list(obsTable)
BM = (product_list['productSubGroupDescription'] == 'IMA')
product_list = product_list[BM]
product_list.show_in_notebook(display_length=5)
Choose a single exposure file to work on. In this example, we choose the first exposure. To create usable data, you will have to follow this work flow on all individual IMA files in your dataset.
myID = product_list['obsID'][0:1]
Download the IMA and FLT files for that exposure. The standard pipeline-FLT will be used for comparison with the detrended final product.
download = Observations.download_products(myID,mrp_only=False,productSubGroupDescription=['IMA','FLT'])
Display the results of the download operation.
download
Read the files that were just downloaded locally. In addition, have the path be just the rootname, i.e. without the file extension.
localpathtofile = download['Local Path'][1][:-8]
localpathtofile
original_ima = fits.open(localpathtofile+'ima.fits')
original_flt = fits.open(localpathtofile+'flt.fits')
original_ima.info()
Print the number of samples and plot the individual reads of the IMA file.
Note: the individual 'SCI' extensions are stored in reverse order, with 'SCI', 1 corresponding to the last read.
nsamp = original_ima[0].header['NSAMP']
print('NSAMP',nsamp)
fig,axarr = plt.subplots((nsamp+3)//4,4, figsize=(9,3*((nsamp+3)//4)))
for i in range(1,4*((nsamp+3)//4)+1):
row = (i-1)//4
col = (i-1)%4
if (i <= nsamp):
immed = np.nanmedian(original_ima['SCI',i].data)
stdev = np.nanstd(original_ima['SCI',i].data)
axarr[row,col].imshow(original_ima['SCI',i].data,clim=[immed-.3*stdev,immed+.5*stdev],cmap='Greys',origin='lower')
axarr[row,col].set_title('SCI '+str(i))
axarr[row,col].set_xticks([])
axarr[row,col].set_yticks([])
else:
fig.delaxes(axarr[row,col])
fig.tight_layout()
Before running reduce_dash, we need to set some environment variables for several subsequent calibration tasks.
We will point to a subdirectory called crds_cache/ using the IREF environment variable. The IREF variable is used for WFC3 reference files. Other instruments use other variables, e.g., JREF for ACS.
os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = './crds_cache'
os.environ['iref'] = 'iref/'
if not os.path.exists('iref'):
os.mkdir('iref')
The code block below will query CRDS for the best reference files currently available for these datasets and update the header keywords to point to these new files. We will use the Python package os to run terminal commands. In the terminal, the line would be:
crds bestrefs --files [filename] --sync-references=1 --update-bestrefs
...where 'filename' is the name of your fits file.
ima_files = glob('*_ima.fits')
for file in ima_files:
command_line_input = 'crds bestrefs --files {:} --sync-references=1 --update-bestrefs'.format(file)
os.system(command_line_input)
Run the DASH pipeline for a single exposure. This procedure showcases the capabilities and customization options of the DASH pipeline.
Note: the following will only work if you are using the notebooks inside of the Notebook directory. wfc3_dash submodule will be properly packaged and installed within the wfc3_tools module sometime in the future.
If you move the notebooks and want to use them elsewhere, you can still provide a temp_path to the dash codes and remove the comments below.
#import sys
#module_path = os.path.abspath(os.path.join('..'))
#if module_path not in sys.path:
# sys.path.append(module_path)
#tmp_path = ".../wfc3_dash/wfc3_dash"
#if tmp_path not in sys.path:
# sys.path.append(tmp_path)
We use the both IMA and FLT extensions of our local image to create a DashData object.
myDash = DashData(localpathtofile+'ima.fits', flt_file_name=localpathtofile+'flt.fits')
print(myDash.root)
A difference (diff) file contains the counts accumulated between two reads of the IMA file. The diff files are written to disk in a directory named ./diff under the current working directory (cwd). In creating diff files from the readouts of the IMA, the first difference, between the 1-st and 0-th read is ignored because of its very short exposure time of 2.9 seconds, resulting in a noisy image.
In order to create the best possible results, the split_ima() method uses the bestrefs function from CRDS to ensure all reference files are up to date and available.
myDash.split_ima()
Print the number of diff files and plot the diff files.
ndiff = len(myDash.diff_files_list)
print('Number of diff files',ndiff)
if ndiff > 4:
fig,axarr = plt.subplots((ndiff+3)//4,4, figsize=(9,3*((ndiff+3)//4)))
for i in range(4*((ndiff+3)//4)):
row = (i)//4
col = (i)%4
if (i < ndiff):
diff_i = fits.open(myDash.diff_files_list[i]+'_diff.fits')
immed = np.nanmedian(diff_i['SCI'].data)
stdev = np.nanstd(diff_i['SCI'].data)
axarr[row,col].imshow(diff_i['SCI'].data,clim=[immed-.3*stdev,immed+.5*stdev],cmap='Greys',origin='lower')
axarr[row,col].set_title('Diff:'+str(i+1))
axarr[row,col].set_xticks([])
axarr[row,col].set_yticks([])
else:
fig.delaxes(axarr[row,col])
else:
fig,axarr = plt.subplots(1,ndiff,figsize=(15,15))
for i in range(ndiff):
immed = np.nanmedian(diff_i['SCI'].data)
stdev = np.nanstd(diff_i['SCI'].data)
diff_i = fits.open(myDash.diff_files_list[i]+'_diff.fits')
axarr[i].imshow(diff_i['SCI'].data,clim=[immed-.3*stdev,immed+.5*stdev],cmap='Greys',origin='lower')
axarr[i].set_title('Diff:'+str(i+1))
axarr[i].set_xticks([])
axarr[i].set_yticks([])
fig.tight_layout()
This file mimics a typical association file for dithered exposures, which is used by AstroDrizzle to align and stack multiple exposures taken at the same sky position with small dithers.
We exploit the fact that a WFC3/IR exposure taken under gyro control can be effectively split into individual pseudo-exposures (the diff images we created in Section 3.2). From there, AstroDrizzle can treat such pseudo-expsoures as individual dithers, and combine them into a single exposure.
myDash.make_pointing_asn()
Show the content of the asn file, which was created in ./diff.
asn_filename = 'diff/{}_asn.fits'.format(myDash.root)
asn_table = Table(fits.getdata(asn_filename, ext=1))
asn_table.show_in_notebook()
Create segmentation map from original FLT image to assist with background subtraction and fixing of cosmic ray flags using create_seg_map(). This method makes a directory called ./segmentation_maps, which holds the outputs.
myDash.create_seg_map()
Plot segmentation map.
rootname = myDash.root
segmap_name = ('segmentation_maps/'+ rootname + '_seg.fits')
segmap = fits.getdata(segmap_name)
fig = plt.figure(figsize=(6, 8))
plt.title(segmap_name)
plt.imshow(segmap, origin='lower', vmin=0, vmax=1, cmap='Greys_r')
Print and read source list.
sourcelist_name = ('segmentation_maps/' + rootname + '_source_list.dat')
sourcelist = ascii.read(sourcelist_name)
print(sourcelist)
Let's create a segmentation map and source list from the difference files. We need to make source lists from our difference files created from the IMA so that TweakReg can better align these difference files to catalogs, each other, etc.
First, generate a list of difference files that contain the full path name.
diffpath = os.path.dirname(os.path.abspath('diff/{}_*_diff.fits'.format(rootname)))
cat_images = sorted([os.path.basename(x) for x in glob('diff/{}_*_diff.fits'.format(rootname))])
sc_diff_files = [diffpath + '/' + s for s in cat_images]
Then, create a difference segmentation map using diff_seg_map() and the diff files.
myDash.diff_seg_map(cat_images=sc_diff_files)
Plot the segmentation map from a diffrence files.
segmap_name = ('segmentation_maps/' + rootname + '_01_diff_seg.fits')
segmap = fits.getdata(segmap_name)
fig = plt.figure(figsize=(6, 8))
plt.title(segmap_name)
plt.imshow(segmap, origin='lower', vmin=0, vmax=1, cmap='Greys_r')
Subtract background from the individual reads taken from the original IMA file using the DRZ and SEG images produced in the background subtraction of the original FLT.
By default, subtract_background_reads() will subtract the background and write it to the header. By setting subtract=False, the background will not be subtracted and will only be written to the header. In addition, setting reset_stars_dq=True will reset cosmic rays within objects to 0 since the centers of the stars are flagged.
myDash.subtract_background_reads()
Now, we can use fix_cosmic_rays() to reset cosmic rays within the segmentation maps of objects and use L.A.Cosmic to find them again.
myDash.fix_cosmic_rays()
Align reads from the IMA to one another by aligning each difference file to the first diff file.
Listed below are all the parameters available to myDash.align(). align() uses TweakReg to update the WCS information in the headers of the diff files, then drizzles the images together using AstroDrizzle. There are more parameters available to users when working with TweakReg and AstroDrizzle that could be an integral part of the workflow for users of DASH. The example below lists the default values set for every input:
myDash.align(self, subtract_background = True,
align_method = None,
ref_catalog = None,
create_diff_source_lists = True,
updatehdr = True,
updatewcs = True,
wcsname = 'DASH',
threshold = 50.,
cw = 3.5,
searchrad = 20.,
astrodriz = True,
cat_file = 'catalogs/diff_catfile.cat',
drz_output = None,
move_files = False)
Refer to documentation to customize parameters for TweakReg and AstroDrizzle.
Note: the error UnboundLocalError: local variable 'sig' referenced before assignment can be solved by lowering threshold parameter.
myDash.align(updatehdr=False, updatewcs=True, astrodriz=False)
Print the shifts file to analyze how well the alignment went. Do not update header until shifts, as seen in the xrms and yrms columns, are satisfactory. Further information about the outputs in the shift file and what is 'satisfactory' can be found in the Drizzlepac Handbook.
shift_file = glob('shifts/shifts_*.txt')
shift_file_name = shift_file[0]
shift_tab = Table.read(shift_file_name,
format='ascii.no_header',
names=['file','dx','dy','rot','scale','xrms','yrms'])
formats = ['.2f', '.2f', '.3f', '.5f', '.2f', '.2f']
for i, col in enumerate(shift_tab.colnames[1:]):
shift_tab[col].format = formats[i]
shift_tab
Let's align our images with a threshold of 20, and update the headers and WCS information.
myDash.align(threshold = 20.)
Plot the final DRZ image and compare to the original IMA.
sci_name = myDash.root + '_drz_sci.fits'
og_flt_name = 'mastDownload/HST/' + myDash.root + '/' + myDash.root + '_ima.fits'
sci = fits.getdata(sci_name)
og_flt = fits.getdata(og_flt_name)
fig = plt.figure(figsize=(9, 4))
ax1 = fig.add_subplot(1,2,2)
ax2 = fig.add_subplot(1,2,1)
ax1.set_title('DASH Pipeline Reduced Science File')
ax2.set_title('Original IMA (not reduced using pipeline)')
ax1.set_xlim(-10,1120)
ax2.set_xlim(-10,1120)
ax1.set_ylim(-10,1050)
ax2.set_ylim(-10,1050)
ax1.imshow(sci, vmin=0, vmax=40, cmap='Greys_r', origin='lower')
ax2.imshow(og_flt, vmin=0, vmax=40, cmap='Greys_r', origin='lower')
plt.show()
Thank you for walking through this notebook. Now using WFC3 data, you should be more familiar with:
wfc3_dash.Below are some additional resources that may be helpful. Please send any questions through the HST Helpdesk.
Authors: Catherine Martlin; WFC3 Instrument Team
Updated on: 2021-10-07
If you use numpy, matplotlib, astropy, astroquery, or drizzlepac for published research, please cite the
authors. Follow these links for more information about citing the libraries below: